Time-Dependent Density Functional Theory for Open Quantum Systems with Unitary Propagation 
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We extend the Runge-Gross theorem for a very general class of Markovian and non-Markovian open quantum 
systems under weak assumptions about the nature of the bath and its coupling to the system. We show that for 
Kohn-Sham (KS) Time-Dependent Density Functional Theory, it is possible to rigorously include the effects of 
the environment within a bath functional in the KS potential, thus placing the interactions between the particles 
of the system and the coupling to the environment on the same footing. A Markovian bath functional inspired 
by the theory of nonlinear Schrodinger equations is suggested, which can be readily implemented in currently 
existing real-time codes. Finally, calculations on a helium model system are presented. 
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Current advances in the manipulation and control of 
nanoscale systems allow for an unprecedented opportunity 
to probe the non-equilibrium dynamics of a wide variety of 
condensed matter systems on a broad range of timescales 
|Ul litl- Serious effort is therefore required for the develop- 
ment of tractable theoretical methods that can shed some light 
on many-body dynamics without directly solving the time- 
dependent Schrodinger equation (TDSE) for an object com- 
posed of many particles. One of the most promising methods 
in this regard is Time-Dependent Density Functional Theory 
(TD-DFT) HUH], which is formally equivalent to the TDSE, 
but is based on the particle density rather than the wavefunc- 
tion. 

Recently, there has been a considerable interest in devel- 
oping an Open Quantum Systems (OQS) formalism for TD- 
DFT, where the number of particles in the system remains 
fixed, but there is energy exchange with an environment 
0, S H H OH 0. This effort allows for the descrip- 
tion of particle transfer within the system, spontaneous de- 
cay, inelastic scattering, and many other ubiquitous relaxation 
and dephasing phenomena. For a Markovian equation of the 
Lindblad form, Burke, Car, and Gebauer (BCG) proved that a 
statement analogous to the Runge-Gross (RG) theorem holds, 
namely, that there is a one-to-one correspondence between the 
time-dependent particle density and the external scalar poten- 
tial provided that the particle-particle interaction, initial quan- 
tum state, and the bath jump-operators remain fixed J3l. To 
place their result in a practical context, BCG assumed the exis- 
tence of a Kohn-Sham (KS) scheme in order to carry out their 
calculations. In the mentioned procedure, an artificial non- 
interacting open system, the so-called KS system, evolves un- 
der an effective KS potential and is expected to reproduce the 
particle density of the original system 01411 . By virtue of their 
theorem, any observable is a functional of the particle density, 
so in principle, the KS system contains all the information 
about the observables of the original system. The question of 
whether or not such a non-interacting KS system exists is not 
obvious, but clearly crucial for KS theory, and it is known as 
the non-interacting V-representability problem 0151 llol . We 



In the case of closed systems, Van Leeuwen has proved 
that it is in fact possible to reproduce the particle density of 
a many-body interacting system with an effective KS poten- 
tial acting on an auxiliary system with no particle-particle in- 
teractions lll5ll . This KS potential is unique, and in general, 
expected to show a nonlinear and nonlocal funcional depen- 
dence on the history of the particle density [17]. Intuitively, 
we can argue that in the KS system we formally give up the 
linearity of the many-body equation of motion for a nonlinear 
surrogate which, nevertheless, is an effective single-particle 
equation. With this in mind, a natural question to ask is: Just 
as with the particle-particle interactions, can we subsume the 
coupling between the system and the bath into an additional 
nonlinearity of the density in the effective KS potential? In 
the next paragraphs, we report that this is indeed the case. 

Consider an -particle open quantum system described by 
a time-dependent density matrix p(f) which, in the position 
representation, is a function of 6N coordinates and time t . The 
most general equation of motion for an open quantum system 
is a master equation of the form (atomic units used through- 
out) JH, 



(>(t) = -i[ti(t),p(t)]+ f \%r(t,t')p(t')dt'+^(t). (i) 

Jo 



Here, H{t) = 



■Li<jU(r h fj) is the gen- 



note that non-interacting V-representability in the context of 
BCG's formalism has been assumed, but not formally proven. 



erator of the unitary piece of the evolution. In general, H(t) 
is an effective renormalized Hamiltonian of the system due 
to its interaction with the bath, where U (fi,rj) is a symmetric 
pairwise interaction potential and V (r,t ) an external scalar po- 
tential. Finally, J(f(t,t r ) is a memory kernel which describes 
the non-unitary effects of the bath on the evolution of the sys- 
tem, and 3?" (t) is an inhomogeneous term which is present 
only if there are initial correlations between the system and 
the bath. Quite generally, J(f(t ,t') and may be functions 
of V(r, ?), such as in the case of a strong laser field interacting 
with a molecule in condensed phase 11911 . whereas 3^(t) may 
also depend on the initial state p (0). 

Furthermore, for notation, we define the operators that mea- 
sure the particle density as n(r) = 8(r — rf) and the current 
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density as j(7) - 
state a theorem. 



j£i{5(r — ri),$i}. We are now ready to 



Theorem.- Let the original system be described by the den- 
sity matrix p(t) which, starting as p(0), evolves according to 
Eq. (Q]i. Consider an auxiliary system associated with the den- 
sity matrix p'(t) and initial state p'(0), which is governed by 
the equation: 

P '(t) = -i[ft'(t),p'(t)]+ [ 'dt'j(r'(t,t')p'(t')+&'(t), (2) 

Jo 

where the functional forms of ,3tf'(t) and 3T'(t) are given, 



and its Hamiltonian reads as H'(t) = -j^-+V (%,t) 

Y*i<jU' (7i,rj), where I/'(r;,r/) is also given. Under mild con- 
ditions, there exists an external potential V'(7, t) that drives the 
auxiliary system in such a way that the particle densities in the 
original and the auxiliary systems are the same at every point 
in time and space, i.e., (h(7))' t = (h(7)) t . This statement is 
true provided that p'(0) guarantees that (n(7))' t=0 = (n(7)) t= o 
and (n(r))' t=0 = (n(7)) t=0 . 



Proof.- We use similar techniques to the ones employed by 
van Leeuwen 1 131 and Vignale 11611 . The detailed steps of a 
related derivation may be found in l20ll . First, by using Eq. 
(Q]l, we can find the equation of motion for the second deriva- 
tive of the particle density of the original system with respect 
to time: 

{h(f)) t = V ■[(h(7)) t VV{7,t)/m + $(r,t) + 

.#(r,r)//n + #(r,r)] + ( /(r,r). (3) 

Here, VV(r,f) is proportional to the external electric 

field, §(?,t) = -lZa,fi&&{l.l{Vta,{*tll,8(?-%)}}) is 
the divergence of the stress tensor, where a,j5 = x,y,z, 
JF(7,t) is the internal force density caused by the pair- 
wise potential ,#(r,f) = 8(7 - r^Zj^i %U(n - 7j)), 

and Sf(r,f) = Tr{j(7)(&dt' 'JT(t ,t')p(t') + 3T{t))} and 
J(7,t) = ^Tr{h(r) %dt' ,X (t j)p(t') + &(t))} are terms 
which arise due to the coupling to the bath. Similarly, by em- 
ploying Eq. (0, it is possible to derive an equivalent equa- 
tion for (h{f))' t , where the variables in Eq. Q are substituted 
by their primed analogues. If we subtract these two equa- 
tions and eliminate the variable (n(7,t))' with the restriction 
(h(7))' t = (h(r)) t , we obtain an identity with time-dependent 
parameters that can be Taylor expanded about t = 0. Denot- 
ing the Taylor expansion coefficients by Ok = n ~ jt^ \i=Qi 
we collect the terms of order t ', and arrive at the expression: 



-V.(„ (r)V(y/ (?))) = 
-V • (m#/(r) + #/(r) + m^'(r)) + m f[{r) 
+ V ■ [m&i if) + ,0, (r) + m$i {r))-m/i (r) 



-v.(« (?)v(y,(r))) + v. t, n ^(yl-k-Vi-k(f)) 

\k=l j 



(4) 



We now make a claim: If the right hand side of Eq. (01 
contains no coefficients V[(r) for k > /, it can be regarded 
as a recursion relation to construct V/ from the lower or- 
der coefficients Vi{f) for < k < I. This would imply that 
each coefficient can be uniquely solved recursively upon the 
specification of a boundary condition, which we can conve- 
niently set to V/(r) — > as \r\ — > °°, for all I. Finally, the ex- 
plicit construction of V'(r, t) through its Taylor coefficients, 
V'(r,t) — Y,kVl(r)t k , proves the theorem. 

If Jf'(t,t') and &'(t) do not depend explicitly on V'(f,t), 
the claim can be systematically shown 1I20I1 . Otherwise, 

^1(7) can depend at most on J^''f 1 = /r ~ ^j''^ |r=o (the 
integral terms f dt(-) naturally vanish at t = 0) and 2F{ — 

T\ ~ l'= 0, General expressions derived with projection- 
operator methods lfl9ll can be used to formally show that 
Jtf'f' should depend at most on V{_ l (7), which supports our 
claim. This fact can be interpreted in very physical terms: the 
action of the external field V'(7,t) on the system is local in 
time through the unitary piece of the master equation. The ef- 
fects of V'(7,t) on the system leak out to the bath and return 
as memory effects through the memory kernel only at times 
t" strictly later than t . In other words, Jtf ! (t,t) can depend on 
V' (?,/') for t' < t, but should not depend on the instantaneous 
V'(7,t). 

A similar conclusion may not be made for arbitrary 
terms, since at t = 0, the initial correlations between the sys- 
tem and the bath may depend on V(r,0), and 2f( could de- 
pend on V!{7). However, as long as depends at most on 
V'j _j (7), the claim and the theorem will necessarily hold. This 
is the only warning of the proof, and this requirement can be 
checked on a case by case basis, but it is easily guaranteed in 
the case of initial factorizable conditions between the system 
and the bath, or if the inhomogeneity is V-independent, which 
occurs if the external field is weak or if the bath is Markovian. 
□ 

Several important corollaries hold from the theorem. If 
p'(0) = p(0), U'{7„7j) = U(7 h 7j), X'itj) = Jf (t,0, and 
^'( f ) = ^( f ),thenEq. ©reads: -V- (noV(V/_ fc -V;_jt)) = 
V • (iLi n k^{y'i- k - Vl-kj) , which means that V( = V, for all 
/. This allows for an extension of the RG theorem to a large 
class of OQS: For fixed initial state, interparticle potential, 
memory kernel and inhomogeneity, there is a one to one map 
between particle densities and scalar potentials. This state- 
ment allows us to regard the time-dependent particle density 
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as a fundamental variable just as the time-dependent density 
matrix. For Markovian equations of the Lindblad form, this 
reduces to the result proven by GCB. 

The theorem also justifies the KS scheme of BCG and its 
generalization to a wide range of OQS, namely, that it is pos- 
sible to choose an auxiliary open system with no particle- 
particle interactions, U'{fi,rj) = 0, to reproduce the same par- 
ticle density as the original system. However, we want to 
take a different approach on the subject and make the obser- 
vation that the proof also allows us to consider the case where 
U'(fi,fj) = X'{t/) = &'{t) = 0, that is, a KS system that 
evolves unitarily as if it were a driven closed system, but still 
reproduces the particle density of the original open system 
that interacts with the bath and evolves through a non-unitary 
equation of motion. Therefore, we have rigorously justified 
the intuition hinted at the beginning of the letter, that is, the 
possibility to conceive of a KS system where we subsume the 
effects of the bath in an additional term in the KS potential 
lE^ll . In this new KS theory, we shall rewrite the KS potential 
as V' — V + Vh + Vxc + Vb a th, where V is the original exter- 
nal potential, V//(r, f) = / cPr 1 is the Hartree term, V xc 
is a standard approximation to the exchange-correlation (xc) 
term due to the many-body effects within the system, such 
as an adiabatic functional [21], and finally, Vbath is the new 
term due to the bath, which includes additional correlations 
on the particles of the system, and which we expect to be non- 
adiabatic. Finally, we must discuss the feasability of the initial 
conditions for our KS scheme. It is always possible to propose 
a pure state single Slater determinant y/'(0) = det[0,-(ry)] 

which satisfies the restriction {xjf' (0)\h(r)\ \j/'(0)) = (n(r)) f= o 
by employing the Harriman construction B22I1 . By defin- 
ing a new state y/(0) = -^det[^>i(fj)e ,ai ^], the set 
of phases {a,} can be chosen with considerable free- 
dom in order to satisfy Jy(i/(f)|«(r)|y/(f))|r=o = — V • 

(l/I^^PVC-jarg^-^ + ai^)) = (h(f))t^o, in which 
case, we can choose y/(0) as the initial KS wavefunction, or 
equivalently, p'(0) = |i/(0))(</(0)| as the initial KS density 
matrix. Note that this argument is irrespective of the purity of 
the initial state of the original system. 

Model system and suggestion of "bath "functional.- We re- 
fer the reader to Ref. 12011 . which reports a numerical study 
that constructs the KS potential V' for a harmonic oscilla- 
tor model coupled to a heat bath. In this letter, we will be 
concerned with the study of a model system, namely, a 1- 
d helium atom 123l 12411 coupled to a heat bath. We write 
the total system-bath Hamiltonian as Hj = H$ + H$b + Hb- 
Hs = E*=i (^ 2 /2 + V(X,-,f)) +W(Xi -Xz) describes the he- 
with Xj and P, denoting the positions and mo- 
W(X) = e 2 /VX 2 + 1 being a soft- 
Coulomb potential, and V (X) = — 2W(X) the external poten- 
tial, which in this case is only due to the nucleus. Hb + Hsb = 



lium atom 
menta of the electrons 



rX 



corresponds to a har- 



trons. We assume that the bath is an infinite set and its 
distribution of couplings can be approximated by a contin- 

2 

uous Ohmic spectral density, J(co) =Y,j 2m a ^i®] ~~ w ) = 

0(c0)^C0e~ ffl /'° l , where 9(co) is the step function, ^ is the 
intensity of the coupling, and co c is a cutoff frequency for the 
bath modes. From a computational point of view, the dynam- 
ics of the composite system-bath object is intractable. Since 
the emphasis is on the system, and not on the bath, we take an 
OQS approach: For weak coupling %q and large CO c , the Born- 
Markov approximation is justified, and it is straightforward to 
obtain a memoryless master equation of the Lindblad form for 
the system. At zero temperature (T = 0), it reads, 



p(0 



-i[H s ,p]-±(VLp 



pL'L- 



2LpV 



(5) 



where H$ = Hs + (x 2 +y 2 ) is a renormalized Hamiltonian 
due to coupling to the bath. We denote \g) and \e) to be the 
ground and first singlet excited states of H$ respectively, so 
that the jump operators L can be expressed in the form L = 
\g) (e\. L promotes quantum jumps from \e) to \g) . The rate of 
these transitions is captured by Y=27i\(e\jj.\g)\ 2 J((O eg ), where 
ll = Y^=\Xi is the dipole operator. 

We proceed to derive a bath functional which could be used 
in the KS theory for TD-DFT applied to systems interacting 
with a Markovian bath, just like our model system. For a 
single particle, Kostin Il25ll has previously constructed a dissi- 
pative nonlinear Schrtidinger equation, where — H\j/, for 

„2 



which the Hamiltonian in 1-D reads H 



2M 



-V + V bath , with 



the bath potential being given by V bath (X , t ) = | In ( ^|jy 
This equation of motion has the very interesting property 
that at the level of observables, it satisfies the Langevin 
equation at T = 0, i.e., (X) = (P) = -A(P) - (^>>, 
as can easily be checked by direct substitution. The friction 
coefficient A may be obtained from a microscopic derivation 
of the Langevin equation, which in the case of a particle 
bilinearly coupled to an Ohmic bath of strength yields 
A = n^o /2. Furthermore, a quick inspection allows us to 
rewrite Vbath as a functional of the particle density ll30"ll . 



V bath [(h(X%,(j(X%](X,t) = A f dX'^^. (6) 



monic bath with bilinear coupling to the positions of the elec- 



(h(x')) t - 



For more than one particle, this identification is not formally 
possible, but regardless, we shall heuristically assume it as our 
Markovian bath functional (MBF) Bill . Non-Markovian gen- 
eralizations of Eq. (|6]l may be readily conceived starting from 
nonlinear Schrodinger equations which reproduce the gener- 
alized Langevin equation for its observables. Physically, this 
suggestion is very appealing: The dragging force due to the 
MBF is proportional to ||pfjjj> which is the velocity field. The 
coefficient A can be approximated from the spectral density 
and conveniently scaled to reflect the many-body coupling to 
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, where a,- = —/arg (</>,). 



the bath. From the single Slater determinant KS wavefunc- 
tion, XjfKsit) = det[&(Xy,f)], we can express V bat h(X,t) = 
I rz dx r Li\UX'j)\ 2 Va,(x',t) 

In order to gain insight on the system of consideration, we 
performed several calculations for which the results are sum- 
marized in Fig. 1. The initial state of the helium atom was 
taken to be the pure state y(0) = ^(Is) + l e ))- We prop- 
agated the system in real time with three different methods. 
For the first method (black solid curve), we evolved the den- 
sity matrix of helium using the master Eq. ©. We chose a 
spectral density with values |o = 0.01 Ef, and a> c = 10E/,. The 
real space eigenbasis of H$ was obtained with the OCTOPUS 
package 12611 . resulting on an energy gap A eg = 0.85 Ef, and 
a dipole moment (e\jj,\g) — \.\a.u. The choice of parame- 
ters justifies the Markovian conditions for the master equation. 
The expected damped oscillations calculated with this method 
are shown in the solid curve. The second method (red solid 
curve) was performed to calibrate the parameter X in V c n s . 
We evolved the time dependent Schrodinger equation with 
the effective Hamiltonian H$ + Vdis using the Suzuki-Trotter 
split operator method lE7tl . where the many -body dynamics 
was computed exactly via Hs, but the coupling to the bath 
entered through the nonlinear dependence of Vdis on («(?)),. 
We scanned several X parameters and found A = 0.075 Ef, 
to reproduce the curve derived from (A) with high accuracy 
1 3211 . Finally, for the third method (black dotted curve) we 
carried out a TD-DFT KS calculation with exact exchange 
and same dissipation rate y as in B, that is, Vks = \Vh + Vdis, 



The result for this last 



with V H (X,t) = vfdX'- , — . 

method yields poor results with unphysical Rabi-like oscil- 
lations. The latter are caused by the absence of correlations 
caused by particle-particle interactions B28I1 . Nevertheless, the 
oscillations decay on a similar timescale to the other calcula- 
tions, and reach a steady state due to the MBF 

In summary, we have formally extended TD-DFT to a large 
class of OQS, and rigurously showed the possibility of in- 
cluding the effects of the bath on the dynamics of the sys- 
tem within a bath functional. The latter enters into a TD-DFT 
calculation on the same footing as the standard exchange- 
correlation functionals exclusively due to many-body dynam- 
ics. We have suggested Eq. (O as the Markovian bath func- 
tional which can be readily implemented in currently existing 
TD-DFT codes. Future work must address the derivation of 
the friction coefficient X from a more systematic procedure, 
the explicit derivation of bath functionals for more complex 
memory kernels, and finally, the role of fluctuations and finite 
temperatures in this unitary propagation KS formalism. 
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lowship program. This work was carried out under the 
DAPvPA contract FA 9550-08-1-0285. 




FIG. 1 . Evolution of the dipole moment of a helium atom coupled to a heat 
bath. We present three different calculations: The black solid curve 
represents the "exact" calculation using a master equation. The red solid 
curve is the propagation of the exact many-body dynamics of helium plus 
the MBF. Finally, the dotted black is the TD-DFT calculation with exact 
exchange and MBF. The last calculation yields poor results due to the 
absence of correlations in the electron interactions. Details of the 
calculations can be found in the text. 
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